Introduction¶
The goal of this project is to introduce methods for modeling Gaussian Processes (GP) whose underlying function values (priors) are treated as latent variables. For the sake of this example, we will be using PyMC3's Gaussian Process function along with the gp.Latent class. This class includes a prior and conditional method that are used to construct random variables representing the unknown function whose distribution is the GP prior or GP conditional. This GP implementation can be used to implement regression on data that is not normally distributed. The GP function, f(x), can be modeled as:
$$
f(x) \sim GP(m(x), k(x, x'))
$$
Where $ m(x) $ is our mean function and $ k(x, x') $ is our covariance (kernel) function.
Both the mean and covariance functions default to zero, however; the covariance function can be set to any of the following:
| Function |
Description |
| Constant(c) |
Constant valued covariance function. |
| WhiteNoise(sigma) |
White noise covariance function. |
| ExpQuad(input_dim[, ls, ls_inv, active_dims]) |
The Exponentiated Quadratic kernel. |
| RatQuad(input_dim, alpha[, ls, ls_inv, ...]) |
The Rational Quadratic kernel. |
| Matern32(input_dim[, ls, ls_inv, active_dims]) |
The Matern kernel with nu = 3/2. |
| Matern52(input_dim[, ls, ls_inv, active_dims]) |
The Matern kernel with nu = 5/2. |
| Exponential(input_dim[, ls, ls_inv, active_dims]) |
The Exponential kernel. |
| Cosine(input_dim[, ls, ls_inv, active_dims]) |
The Cosine kernel. |
| Periodic(input_dim, period[, ls, ls_inv, ...]) |
The Periodic kernel. |
| Linear(input_dim, c[, active_dims]) |
The Linear kernel. |
| Polynomial(input_dim, c, d, offset[, ...]) |
The Polynomial kernel. |
| WarpedInput(input_dim, cov_func, warp_func) |
Warp the inputs of any kernel using an arbitrary function defined using Theano. |
| Gibbs(input_dim, lengthscale_func[, args, ...]) |
The Gibbs kernel. |
| Kron(factor_list) |
Form a covariance object that is the kronecker product of other covariances. |
| Coregion(input_dim[, W, kappa, B, active_dims]) |
Covariance function for intrinsic/linear coregionalization models. |
Priors¶
The prior method places a multivariate normal prior distribution on the vector of function values, $ f $. This can be written as:
$$
f \sim MvNormal(m_x, K_{xx})
$$
where the vector $ m $ and the matrix $ K_{xx} $ are the mean vector and covariance matrix, respectively, evaluated over the inputs $ x $.
Since the function variables are considered to be latent, we aren't going to be explicitly setting them. Instead, we will initialize the mean and covariance functions. An example of this is shown below:
import numpy as np
import pymc3 as pm
# A one dimensional column vector of inputs.
X = np.linspace(0, 1, 10)[:,None]
with pm.Model() as latent_gp_model:
# Specify the covariance function.
cov_func = pm.gp.cov.ExpQuad(1, ls=0.1)
# Specify the GP. The default mean function is `Zero`.
gp = pm.gp.Latent(cov_func=cov_func)
# Place a GP prior over the function f.
f = gp.prior("f", X=X)
The last line of this code block is evaluating the initialized functions across all of the inputs, X, to return a vector of function values, $ f $. It is important to note that the GP function requires X to a one-dimensional column vector of inputs.
When running the prior across inputs X, PyMC3 also does some additional processing in the background. PyMC3 will reparameterize the prior by rotating it using Cholesky matrix factorization. Cholesky factorization decomopses a matrix into the product of a lower triangle matrix and its conjugate transpose -- which is useful for efficient numerical solutions. The resulting reparameterized model is:
$$
v \sim N(0,1) \\
L = Cholesky(L_{xx}) \\
f = m_x + Lv
$$
What is a "lower traingular matrix"? Assume we have some matrix, L, the lower triangle matrix would be all of the values to the left of, underneath, and including the diagonal. See example below.
$$ L =
\begin{bmatrix}
l_{1,1} & & & & 0 \\
l_{2,1} & l_{2,2} \\
l_{3,1} & l_{3,2} & \ddots \\
\vdots & \vdots & \ddots & \ddots \\
l_{n,1} & l_{n,2} & \dots & l_{n,n-1} & l_{n,n}
\end{bmatrix}
$$
What is the "conjugate transpose"? This is the process of first taking the transpose of a matrix, basically rotating it 90 degrees in a sense, and then taking the conjugate of every entry in the matrix. In our case, we actually need to take the complex conjugate, which introduces imaginary numbers:
- Assume we have the matrix $ A $
- $ A = \begin{bmatrix} 1 & -2-i & 5 \\ 1+i & i & 4-2i \end{bmatrix} $
- We first transpose (rotate) the matrix:
- $ A = \begin{bmatrix} 1 & 1+i \\ -2-i & i \\ 5 & 4-2i \end{bmatrix} $
- Then we conjugate every element:
- $ A = \begin{bmatrix} 1 & 1-i \\ -2+i & -i \\ 5 & 4+2i \end{bmatrix} $
The resulting decomposition is typically written in the form $ A = LL^* $, where $ L $ indicates the lower triangular matrix and $ L^* $ is the conjugate transpose. See this link for more information.
Conditional¶
The conditional method implements the predictive distribution for function values that were not part of the original dataset. This distribution is,
$$
f_* | f \sim MvNormal(m_* + K_{*x}K_{xx}^{-1}f,K_{**}-K_{*x}K_{xx}^{-1}K_{x*})
$$
This distribution is created using the values $ f $ that the GP prior was over, and is now the set of new points $ f_* $. Using this we can can construct a random variable with:
# vector of new X points we want to predict the function at
X_star = np.linspace(0, 2, 100)[:, None]
with latent_gp_model:
f_star = gp.conditional("f_star", X_star)
Key Takeaways¶
- We use
pymc.gp.Latent for unknown function variables
pymc.gp.Latent.prior allows us to fit a model to a given datasetpymc.gp.Latent.conditional allows us to predict and forecast values
